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This paper examines the interplay of the effect of cross immunity and antibody-dependent en- 
hancement (ADE) in mutistrain diseases. Motivated by dengue fever, we study a model for the 
spreading of epidemics in a population with multistrain interactions mediated by both partial tem- 
porary cross immunity and ADE. Although ADE models have previously been observed to cause 
chaotic outbreaks, we show analytically that weak cross immunity has a stabilizing effect on the 
system. That is, the onset of disease fluctuations requires a larger value of ADE with small cross 
immunity than without. Ifowever, strong cross immunity is shown numerically to cause oscillations 
and chaotic outbreaks even for low values of ADE. 
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The spreading of infectious diseases having mul- 
tiple strains in a population can exhibit very com- 
plex dynamics, ranging from periodic and quasi- 
periodic outbreaks to high dimensional chaotic 
behavior. Several sociological and epidemiolog- 
ical factors characterize the disease spread at dif- 
ferent levels, such as interactions among the dis- 
ease strains, social contacts, and human immune 
responses. In this work we focus on dengue fever, 
a vector born disease which has exhibited as many 
as 4 different strains, and is endemic in large ar- 
eas of Southeast Asia, Africa and the Americas. 
A notable feature of dengue is its interaction with 
the human immune system. When an individual 
is infected with dengue, the immune system trig- 
gers an antibody response which will temporarily 
protect against secondary infections. However, 
when the level of protection decreases, secondary 
infections may be possible and the presence of 
low level antibodies triggers an increase in the 
infectiousness of the individual. This effect is 
called antibody-dependent enhancement (ADE). 
In this paper we study a mathematical model for 
the spreading of dengue fever. While ADE alone 
is proved to trigger large amplitude chaotic oscil- 
lations, we show that including weak temporary 
cross immunity stabilizes the system. In contrast, 
we also show that strong cross immunity destabi- 
lizes the dynamics. These results will help under- 
stand implementation of proper control strategies 
when using future vaccines. 
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I. INTRODUCTION 



Understanding the dynamics of muhistrain diseases is 
a key topic in population biology. A suitable model class 
for such diseases, which include influenza, malaria and 
dengue [ll] , must take into account the possibility of in- 
teractions among the serotypes or strains. The nature 
of multistrain interactions strongly affects the impact of 
the disease on the population as well as the mechanisms 
for its control. 

One prominent example of an endemic multistrain dis- 
ease is that of dengue and dengue hemorrhagic fever 
(DHF). Located in Africa, the Americas, and Southeast 
Asia, dengue is one of several emerging tropical diseases 
[31 . There is no vaccine, although clinical trials are un- 
derway in order to generate an immune response across 
all strains Approximately 2.5 billion people are at 
risk for contracting dengue [3l|, (sj , and between 50 and 
100 million cases are reported each year The dom- 
inant four dengue viruses have progressively spread ge- 
ographically to virtually all tropical countries to create 
a global pandemic resulting in several hundred thousand 
hospitalizations every year p^ . Since dengue is so far 
reaching and endemic, it is important to understand how 
it fluctuates in time, so that when proper vaccines are de- 
veloped, implementation may be guided by a more thor- 
ough understanding of the disease. Dengue is known to 
exhibit as many as four coexisting serotypes (strains) in 
a region such as Thailand. Dengue displays a distinc- 
tive mechanism of interaction among the strains, called 
antibody-dependent enhancement. Once a person is in- 
fected and recovers from one serotype, life-long immu- 
nity to that serotype is conferred. Antibodies are de- 
veloped specifically for the first challenging serotype and 
not the other serotypes. In the presence of a new sec- 
ondary infection, low level antibodies developed from the 
first infection form complexes with the second challenging 
serotype so that the virus can enter more cells, increas- 
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ing viral production Viral loads are associated with 
transmissibility, and it is hypothesized that individuals 
with secondary infection are more infectious than dur- 
ing their first infection. This increased transmission rate 
in subsequent infections is known as antibody-dependent 
enhancement (ADE). In vitro studies of dengue fever sug- 
gested that the ADE phenomenon may be due to the in- 
creasing of the infection of cells bearing the IgG receptor 
(G-immunoglobulin) j20| . 

The impact of ADE on the modeling of multistrain dis- 
eases such as dengue is quite profound 12]. In general, 
the first models were of SIR type, with ADE included, 
and they showed that for sufficiently high ADE, oscilla- 
tions were possible. In contrast, single strain SIR models 
only have isolated equilibria and cannot show fluctua- 
tions without external seasonal drives or noise. Recent 
work has begun to analyze in detail the effect of ADE 
quantitatively on the dynamics (H, |2^, as well as the 
competition between serotypes [9(]. It is also still unclear 
if ADE increases transmission of the disease or increases 
mortality, shortening the effective infectious period. The- 
oretical studies suggested that the former case allows for 
coexistence of strains with periodic and chaotic disease 
outbreaks while in the latter the phenomenon may 
decrease persistence [1^ . Throughout this work we shall 
assume the first case to hold. 

In addition to ADE, another type of interaction be- 
tween the strains occurs. Recently, cross protection, or 
cross immunity between serotypes, has been conjectured 
to play a role in the dynamics of dengue [1, . While 
a primary dengue infection with a particular serotype 
may confer long-life immunity to that strain [2l|, |23| , it 
may also confer temporary cross immunity to the other 
serotypes. Cross immunity may act like a prophylactic 
to different strains and may also possess different effica- 
cies. That is, cross immunity may be total (i.e., a period 
of cross immunity always occurs in the event of a pri- 
mary infection), or partial (i.e., only a fraction of the 
actual infected population becomes cross immune before 
being susceptible again to the disease). In general, the 
length of the cross immunity period may vary depend- 
ing on the disease. Cross immunity may result from an 
immunological response to the disease. It acts to reduce 
the susceptibility to a secondary infection, lowering the 
effective probability for reinfection to happen In the 
case of dengue fever, cross immunity may last from two 
up to nine months ^30i] , after which the antibodies have 
dropped to sufficiently low levels that allow infection with 
other strains and subsequent ADE. Cross immun ity p lays 
a crucial role in the co-circulation of strains [1, l26j | and 
the pathogen diversity [l|, . 

Several studies involving separately cross immunity 
P, i, M, \M and ADE j., 1, 12, 28.] have been published 
in the past. The presence of both ADE and cross im- 
munity in such models has not been extensively studied, 
although some recent models have begun to address this 
interaction. As an example of such a model, Ref. stud- 
ied the impact of ADE on the dynamics of a multistrain 



disease with temporary cross immunity, giving partic- 
ular importance to the "inverse ADE" hypothesis (i.e., 
reduced infectivity of secondary infections). Ref. [1] con- 
sidered different types of ADE while allowing for lifelong 
partial cross immunity. Ref. 30] showed that including 
both ADE and temporary cross immunity is necessary 
to produce periodicities consistent with epidemiological 
data. Finally, Ref. [24] included different mechanisms of 
cross immunity in a model with ADE in order to test the 
impact of a period of cross protection on the incidence 
of secondary dengue cases. They found that including 
clinical cross immunity, in which a challenge with a pre- 
viously unexperienced serotype results in an increase of 
immunity towards the challenging serotype, gives inci- 
dence patterns of secondary dengue infections that are 
compatible with collected data. 

The aim of our work is to study in detail the impact 
of both ADE and temporary, partial cross immunity on 
the dynamics of the mutistrain diseases. The outline of 
the paper is as follows: In Section |TT] we introduce the 
model, in Section Hill we analyze the effect of weak cross 
immunity on the system, in Section IIVI we restrict our- 
selves to the case of no ADE to investigate the impact of 
strong cross immunity on the dynamics, and in Section 
|V]we include also ADE and study the interplay between 
cross immunity and ADE. Section IVII concludes with a 
summary and discussion. 



II. DESCRIPTION OF THE MODEL 



The dynamical system considered in this paper is based 
on the SIR (Susceptible - Infected - Recovered) model 
and is a generalization of a multistrain model with ADE 
studied previously 0, HI]. We write the model for an 
arbitrary number n of serotypes, and we include both 
ADE and cross immunity in the dynamics. A set of or- 
dinary differential equations describes the rate of change 
of the population in each of the classes. We assume the 
population size to be normalized to unity, so each state 
represents a fraction of the total population. The quanti- 
ties that enter the equations are the fraction of suscepti- 
bles to all serotypes, denoted by s; the primary infectives 
with strain i, xf, the cross immunes that are recovered 
from strain i and have temporary cross immunity to all 
strains, cf, the recovereds from strain i that are immune 
to strain i only, rf, and the secondary infectives with 
strain j previously infected with strain i =^ j, Xij. The 
flow of an individual through the population in the two 
strain case is shown in Figure [1] Tertiary infections are 
not included 0], so all individuals enter the completely 
immune class rtot after recovery from a secondary infec- 
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tion. The dynamical system is as follows: 



ds 
dt 

dxi 
'dt 

dcj 
'dt 



dn 
dt 

dx ij 
"df 




TABLE I: Parameters used in the model 
Parameter Value Reference 



jjL, 1/host lifespan, years"^ 


0.02 




/3, transmission coefficient, years" ^ 


200 


[12] 


a, recovery rate, years"^ 


100 


[15] 


9, rate to leave the cross 


2 


[30] 


immunity compartment, years'^ 






(j), ADE factor 


> 1 




e, strength of cross immunity 


- 1 





-axi 



l^dXij 



where the parameters are the number of strains n, the 
contact rate (3, the recovery rate cr, the ADE factor 0, 
the strength of cross immunity e, the rate for cross im- 
munity to wear off 9, the birth rate /x, and the mortality 
rate iid- The model of Eqs. [1] allows for one reinfec- 
tion. The parameter e determines how susceptible the 
cross immune compartments Ci are to a secondary infec- 
tion, where e = means no cross immunity effect (the 
Ci compartments are infected as easily as the recovered 
compartments r^) and e = 1 confers complete cross im- 
munity [ci are not susceptible to a secondary infection). 

I3{1 - e) 



(1) 



rate to be either /i^ = /i to maintain a constant popula- 
tion or /id = in our analytical approximation for ease 
of analysis. Parameter values are summarized in Table[Tl 
We vary the ADE (j) and cross immunity strength e as 
bifurcation parameters. 

\ The case without cross immunity, e = 0, reduces to a 
Xkj Previously studied model with only ADE 0, HI] because 
/the cross immune and recovered compartments have the 
same infection rate and are treated identically. It has 
been shown 0, HI] that as ADE is increased, the system 
undergoes a Hopf bifurcation to stable periodic oscilla- 
tions and then to chaos (Fig. ^ . Desynchronization be- 
tween strains occurs in the regions of chaotic outbreaks, 
but all strains are synchronized near the Hopf bifurcation 
when the outbreaks are periodic. The system has been 
analyzed in the neighborhood of the Hopf bifurcation us- 
ing a reduced model that assumes a lower dimensional, 
synchronized system. In the next section, we extend this 
analysis to the case of weak cross immunity. 




rtot 



I3{1 ~ ej 

FIG. 1: Flow diagram of how an individual would proceed 
through the model in the case of 2 serotypes. Note the re- 
duction of susceptibility to a secondary infection through the 
cross immunity factor (1 — e) and the enhancement of sec- 
ondary infectiousness due to the ADE factor cj>. Death terms 
for each compartment are not included in the graph for ease 
of reading. 



Throughout this paper, 

-1 



use 

-1 



= 4 serotypes, 

(3 = 200 years^^, a = 100 years^^, 9 — 2 years""'^, and 
fi = 0.02 years^^ in all numerical simulations [1^. The 
parameter 9~^ is the average time span of cross immu- 
nity, which typically ranges from 2 to 9 months [s^l . We 
choose 9 = 2 years" ^, corresponding to 6 months of cross 
immunity, but we have used 9 = 4 years" ^, equivalent to 
3 months of cross immunity, with no significant difference 
in the results. For convenience, we choose the mortality 




FIG. 2: Bifurcation diagram in ADE for the multistrain sys- 
tem with no cross immunity. For each ADE value, we show 
the local maxima (black) and local minima (gray) of the sus- 
ceptibles during a 100 year time series, after removal of tran- 
sients. From 12811. 



III. STABILIZING EFFECT OF WEAK CROSS 
IMMUNITY 

We consider first the effect of weak cross immunity, 
e ^ 1, and show that it helps stabilize the steady state. 
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Numerical simulations indicate that when e is small, the 
system undergoes a Hopf bifurcation as ADE is increased, 
as it does for the system without cross immunity, and 
the compartments are identical across all n strains near 
the Hopf bifurcation. Thus the system's dimensionality 
reduces. Assuming symmetry between strains, we rewrite 
Eqs. [T] as follows: 



using Eq. [5l four of the five eigenvalues can be obtained: 



dyi 
dt 

dy2 
dt 

dm 
dt 

dm 
dt 

dy5 
dt 



/i - /3nyij/2 - (3(j)n{n - l)yiy^ 
Pvm + P(l){n - l)yiy5 - cry2 

<jy2 - f3{l - e){n - I)y2y3 (2) 
-/3(1 - 6)<^(n - I)2y3y5 - eys 

%3 - P{n - I)y2y4 - f3(t){n - I)^y4y5 

- e)y2y3 + - e)0("- - 1)2/32/5 + /32/22/4 
+/3(/)(n - 1)2/42/5 - crj/5 



where 2/1 represents the fraction of the population that 
is susceptible to the disease, y2 the primary infectives, 
2/3 the cross immunes, j/4 the recovereds, and 2/5 the sec- 
ondary infectives. Since fj,d is a small parameter, for ease 
of analysis we set the mortality rate /id = 0. This ap- 
proximation is equivalent to assuming that all mortality 
occurs in the rtot class, those who have recovered from 
infections with two serotypes. In a region where dengue 
is very common and dengue infections occur early com- 
pared to the human life expectancy, it may be an accurate 
assumption. The endemic steady state for Eqs. [5] is 



2/1 

2/2 
2/3 



(l + 0)/3 
a n 



(3n{n - - e){l + (jj) + 0an 



(3) 



2/4 



2/5 



f3 {l + (t>) [/3/i(l + 0)(1 - e)(n - 1)2 + a0n{n - 1)] 



an (71 — 1) . 



Evaluating the eigenvalues of the Jacobian of Eqs. [2] 
at the steady state allows us to study its stability as 
a function of e. Since both /i and e are small parame- 
ters, we expand the root of the characteristic polynomial, 
P{x{p, e)), of the Jacobian matrix as follows: 

x{fi, e) = Xf) + xifi + X2e + x^fi^ + x^e^ + x^^e. (4) 

Let us also use the following transformation for the char- 
acteristic polynomial P{x) 

P{x) = ^iP{x). (5) 
Substituting Eq.|3]into the characteristic polynomial and 



Ai 
A2 

^3/4 



an 



1 



- 1) 
2^(6*2-^/3)' 



(6) 
(7) 

(8) 



^[02(„_l)_n(0+l)]/i 



pe(l){n - 1) 
2n(e'2 + 13) ' 



The last eigenvalue can be obtained by performing the 
following substitution in the characteristic polynomial 



P'{x) = tj'Pix/ii). 
The fifth eigenvalue is then found to be 
Ar ~ -a. 



(9) 



(10) 



The real part of the pair of complex eigenvalues A3/4 
determines the stability of the system, since the other 
eigenvalues are clearly negative. Notice that the param- 
eter e occurs in both the real and imaginary parts of the 
eigenvalues. Therefore, we expect that e will modify not 
only the stability of the endemic state but also the ensu- 
ing frequency of oscillations. To first order, the real part 



of A 



3/4 



K[A3/4] 



[(t)'^{n-l)~n{(j) + l)\ 



ei34>{n - 1) 

M02 



2an ^' '"'"^ ' ^'^ 2n(6l2-h/3) 

(11) 

Notice that the onset of a Hopf bifurcation is clearly a 
function of /i, e, and (p. By visual inspection of Eq. 111! 
we see that when e is increased from 0, the eigenvalue 
becomes more negative, so cross immunity is stabilizing 
in the limit of small e and /i. 

In Figure [3l we plot the zeros of Equation [Tl] in 0-e 
space, showing the predicted location of the Hopf bifur- 
cation in the presence of ADE and weak cross immunity. 
Below the curve, the steady state is stable. As the cross 
immunity is increased, a larger ADE value is needed to 
destabilize the steady state. Thus weak cross immunity 
is stabilizing. Figure [3] also shows the actual location of 



These were computed 
Note that the Hopf 



the Hopf bifurcation for Eqs. 
using a continuation routine 
bifurcation in Fig. [31 where both the numerical and the 
analytical curves were obtained in the case of no mortal- 
ity, occurs at a larger value of ADE than in the system 
with mortality. However, the predicted trend of stabi- 
lization due to cross immunity is observed in cither case. 
(C.f. Figures 9-10.) 



IV. CROSS IMMUNITY AS CRITICAL 
PARAMETER 

We next study numerically the effect of stronger cross 
immunity. We first consider the case of no ADE (</> = 1) 
and fix the number of strains to n = 4 as for dengue. We 
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Unstable steady state 
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Stable steady state 
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e (Cross Immunity) 





FIG. 3: Predicted and actual location for Hopf bifurcation 
as a function of e and for weak cross immunity in the case 
of no mortality (pd ~ 0). The full curve is the analytical 
prediction (zeros of Eq. Ilip . while the dashed curve is the 
actual location of the Hopf bifurcation, obtained numerically 
for the full system in the case of no mortality. The number 
of strains is n = 4, and other parameters are as listed in the 
text. 

introduce partial cross immunity by increasing the value 
of e continuously from e = (no cross immunity) to e = 1 
(complete cross immunity). The attracting bifurcation 
structure is depicted in Fig. ID 

0.36 I , , , , 1 

0.35 - 



0.34 - 




0.27 I ' ' ' ' ' 

0.2 0.4 0.6 0.8 1 

£ 

FIG. 4: Bifurcation diagram for the system of ODEs (Eq. [!}, 
in absence of ADE {4> ~ 1). The cross immunity parameter e 
is varied from to 1. For each cross immunity value we plot 
the maxima (black) and the minima (gray) of the susceptibles 
during a 100 year time series, after removal of a transient. A 
transition to chaos occurs at e « 0.2. 

For weak cross immunity, the endemic steady state is 
stable. A loss of stability occurs at eH — 0.165. Numer- 
ical analysis of the eigenvalues of the Jacobian of Eqs. [1] 
at the steady state shows that a super-critical Hopf bi- 
furcation occurs, and simulations show that the periodic 
orbit that appears just past the Hopf point is stable over 
a very small range of e. The strains are desynchronized 
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FIG. 5: Poincare section showing X2{n + 1) vs X2{n). e = 

0. 179, = 0. See text for details. 

on the periodic branch, so it is not possible to analyze 
this bifurcation using a reduced model as in Section Um 
Also in contrast to the Hopf bifurcation for weak cross 
immunity studied in the previous section, for which one 
complex pair of eigenvalues loses stability, at the Hopf bi- 
furcation eH three identical complex pairs of eigenvalues 
become unstable simultaneously [13]. 

For en < e < Cc, where ec « 0.20, the system dis- 
plays quasiperiodicity. FigurelHlshows a Poincare map for 
e = 0.179, where the system is quasiperiodic. The map 
is obtained as follows: in the n— dimensional phase space 
an 71 — 1— dimensional surface is introduced by fixing 
the value of one of the variables, in this case the num- 
ber of primary infectives currently infected with strain 

1, xi. We then sample the other variables every time 
their path crosses the hyper-plane, that is, every time xi 
is identical to a fixed value. If the system is periodic, 
then the Poincare map would result in a point, whereas 
if the system is quasiperiodic we obtain a closed curve, 
which is indeed what happens for the times series of Fig- 
ure [5] We have observed two attracting quasiperiodic 
attractors with overlapping regions of stability. Sample 
time series for the quasiperiodic attractors are shown in 
Figure [ni^a),(c). The four strains are desynchronized on 
the quasiperiodic attractors, but with different phase dy- 
namics. 

To study the strain desynchronization in more detail, 
we define a phase difference between compartments, as 
in [2I]. Let Y{t) be the reference compartment and 
Z{t) another compartment. Let {tk} denote the se- 
quence of times for local maxima of Y{t) and {rfc} 
the sequence of times for local maxima of Z{t). For 
Tm G {tk,tk+i), define the phase of Z relative to Y as 
"^ZYiTm) — 2-K ^'^~^'f . The phases of the other primary 
infective compartments relative to xi for the quasiperi- 
odic attractors are shown in Figure [ni[b),(d). For the 
attractor at weaker cross immunity, the phases of the 
strains relative to each other are approximately constant. 
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FIG. 6; Quasiperiodic attractors for e = 0.179 (panels a,b) 
and e = 0.2 (panels (c),(d)). Time series of primary infectives 
(log variables) are shown in (a),(c). Phase differences of pri- 
mary infectives X2,xz,X4, relative to primary infective x\ are 
shown in (b),(d). The time series in (a),(c) are the beginning 
of those used to generate (b),(d). The reference strain x\ is 
the lightest gray curve in (a),(c). Other parameters: = 1. 



This is sometimes called a splay phase state in the cou- 
pled oscillator literature. In contrast, the behavior at 
stronger cross immunity is more complex and qualita- 
tively different, with the order of the strain outbreaks 
changing over time. Finally, since all the strains have 
identical parameters, we note that any permutation of 
strain labels gives another similar quasiperiodic state. 

When the cross immunity is increased above Cc — 
0.20, the system bifurcates to chaos. The presence 
of chaos in SIR multistrain models with cross immu- 
nity has been already revealed by several studies in the 
past [l,[Ill[ll,[23|. We have confirmed the chaotic behav- 
ior by computing the maximum Lyapunov exponent for 
Eqs.[T] The maximum Lyapunov exponent was obtained 
by integrating the linear variational equations along solu- 
tions to Eqs. [T]for 10* years after removal of transients. 
Results are shown in Fig. [71 For e < en, the endemic 
steady state is stable and the maximum Lyapunov ex- 
ponent is negative. For e e (£h,£c), the system exhibits 
quasiperiodic solutions and the maximum Lyapunov ex- 
ponent is zero. For e > ec, the system is chaotic and 
positive Lyapunov exponents are observed. 

Sample time series for chaotic solutions are shown in 
Fig. [5Ja),(c). Panel (a) shows the four primary infective 
compartments, which are desynchronized. We measured 
the phase differences of the other primary infectives rel- 
ative to primary infective xi, and they are frequently 
nonzero, although there appears to be some structure 
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FIG. 7: Maximum Lyapunov exponent of Eqs. [T] for = 1 
as a function of cross immunity strength e. Equations were 
integrated for 10* years after removal of transients. 



with certain phase differences more probable than oth- 
ers, as shown in Fig. [SUb). Figure [DJc) shows times se- 
ries of all primary and secondary infective compartments 
that are currently infected with strain 1. We observe 
that primary and secondary infective compartments in- 
fected with the same strain (i.e., Xi and the three Xj^i 
compartments, where j ^ i) are usually synchronized. 
Figure [Sfd), a histogram of phase differences of the xji 
relative to xi, shows the synchronization more clearly. 
This effect has been observed previously for the model 
with ADE only ^] and has been explained by a collapse 
of the dynamics onto a lower dimensional center mani- 
fold [29^. The same reduction in dimension is observed 
in the system with cross immunity (Figure [8]) as well as 
in the system with both ADE and cross immunity (data 
not shown). 



V. INTERACTION OF STRONG CROSS 
IMMUNITY AND ADE 

We now turn to the interaction of both ADE and cross 
immunity, computing bifurcation diagrams using a con- 
tinuation routine [lO| . Figure O shows the full bifurca- 
tion diagram in (jj-e space. Here, the cross immunity 
ranges from to 1. The vertical axis is a logarithmic 
scale for (f>. The curves show the parameters of (e, 0) at 
which a Hopf bifurcation occurs. However, the curves 
denote different types of stability exchange. When cross- 
ing the black curve, only one pair of eigenvalues crosses 
the imaginary axis, indicating a simple bifurcation to or 
from periodic orbits. In contrast, when crossing the gray 
curve, the situation is degenerate in that 3 identical pairs 
of eigenvalues cross the imaginary axis. In this case, it 
is expected that complicated dynamics such as torus bi- 
furcations may come into existence. For example, when 
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FIG. 8: Chaotic attractor for e = 0.35, 0=1. (a) Time series 
of all four primary infectives (log variables). Black dashed 
curve: xi; darkest gray curve: X2- (b) Histogram of phase 
differences of primary infective X2 relative to primary infec- 
tive xi. (c) Time series of the first primary infective xi and 
the secondary infectives currently infected with strain 1 
(log variables). Black dashed curve: xi. (d) Histogram of 
phase differences of secondary infective X2,i relative to pri- 
mary infective xi. Phase difference histograms are collected 
for 2000 year time series. 



traversing regions i, ii, and vi by increasing (j), we go 
from steady state through a periodic orbit and possibly 
aperiodic behavior, and then through another Hopf bi- 
furcation to return to a steady state. On the other hand, 
if we go from region iii to v by increasing 0, we go from 
periodic or aperiodic behavior through a reverse degen- 
erate Hopf bifurcation to steady state. 

Notice that there are two relatively large regions of 
stable steady behavior: one for small e and small (j) in 
region i, and one for large e and large (f> in region v. (Note 
that the latter region of stable endemic states extends to 
small e and large </), labelled region vi in the figure.) For 
large cross immunity where e is between 0.65 and 1, a 
sufficiently large value of (j) will stabilize the steady state 
state again. However, the value of </> is so large (the 
Hopf bifurcation has values of on the order of 100) 
that it is unrealistic. Therefore, to explore in more detail 
the bifurcations occuring at reasonable values of </>, we 
examine the case where e is small, which is shown in 
Fig. [101 

In Fig. 1101 there are four distinct regions describing 
the stability of the steady state behavior. In region I, 
the endemic steady state is stable. The solid curve is a 
line of Hopf bifurcations where one complex pair of eigen- 
values becomes unstable. The dashed curve is a line of 
Hopf bifurcations where three identical complex pairs of 
eigenvalues become unstable. Therefore, the system has 
zero unstable eigenvalues in region I, one unstable pair 
in region II, three unstable pairs in region III, and four 



FIG. 9: Full bifurcation diagram in cross immunity e and 
ADE (j>. Curves indicate location of Hopf bifurcations. See 
text for details. 



unstable pairs in region IV. At the Hopf bifurcation be- 
tween regions I and II, a stable periodic orbit emerges 
for which all four strains are synchronized and identi- 
cal. This bifurcation was studied in Section IIIII using 
a reduced model that assumed symmetry between the 
strains. This periodic orbit has a narrow region of stabil- 
ity, and then it quickly bifurcates to chaos, so the major- 
ity of region II displays chaotic dynamics. At the Hopf 
bifurcation between regions I and III, the region of sta- 
ble periodic orbits is even smaller, and then the system 
goes to a quasiperiodic attractor. When e becomes suf- 
ficiently large, the system bifurcates to chaos. Although 
quasiperiodic orbits are observed for portions of region 
III shown in Figure 1101 the majority of region III for 
e > 0.2 displays chaotic dynamics. Chaotic dynamics are 
also observed in most of region IV. Figure [10] (inset) also 
partially explores the sensitivity of the average oscillation 
period in region II with respect to e. Here (/> « 3.877, and 
we vary e to compute a branch of periodic orbits. Plotted 
is the period of the branch of periodic orbits (unstable). 
Notice that in the linear range near the bifurcation point 
where e G (0.05,0.07), the slope is on the order of 100, 
showing a clear sensitive dependence of the oscillation pe- 
riod on the cross immune response. For larger values of 
e, the period exhibits a nonlinear response at the turning 
point, resulting in a bi-unstable branch of periodic orbits. 

Finally, we show the interplay between ADE and cross 
immunity by comparing the bifurcation diagram of Fig- 
ure [H which was obtained with a model equivalent to the 
model of Eqs. [1] with no cross immunity, to the case of 
weak and strong cross immunity. In other words, we fixed 
the value of e and built the bifurcation diagram using 
as critical parameter. 



Figure 



11 a 



shows the effect of 
the inclusion of weak cross immunity (e = 0.05). By vi- 
sual comparison with Figure [H it is clear that the region 
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(a) 




FIG. 10: Blowup of bifurcation diagram in Fig. [9] for small 
cross immunity e and ADE <j!>. Curves indicate location of 
Hopf bifurcations. Only region I has stable steady states. 
The inset shows the period of a branch of unstable orbits as 
a function of e for 4> « 3.877 in region II. See text for details. 



of stability is increased: A value of w 2.5 is needed to 
destabilize the system in comparison with w 1.7 needed 
in the case of no cross immunity. Figure [11(b) shows the 
effect of strong cross immunity on the bifurcation struc- 
ture (e = 0.6). The system is observed to be chaotic for 
all the considered values of ADE. 



VI. CONCLUSIONS AND DISCUSSION 

In this work we analyzed the impact of two types of 
strain interactions in a multistrain model for epidemics, 
cross immunity and ADE. The ADE parameter measured 
an increase in infectiousness of secondary infectives, and 
the cross immunity strength determined the reduction in 
susceptibility to other strains during a temporary period 
after recovering from primary infection with one strain. 

The nature of the observed dynamics depended on the 
strength of the cross immunity. Weak cross immunity 
was found to stabilize the endemic steady state. This 
effect was motivated analytically by studying a reduced 
model for weak cross immunity with symmetry between 
strains. Although the analysis was performed for a per- 
turbed system without mortality, both the analytical 
treatment and numerical simulations of the full system 
were in good qualitative agreement. Since the onset of 
fluctuations is determined by Hopf bifurcations in models 
for dengue, the stabilizing effect of cross immunity shows 
that it is an important parameter to include when mod- 
eling disease fluctuations about equilibria. In addition, 
since cross immunity has a strong effect on the period of 
oscillation, it will play a role in determining the timing 




1 1.2 1.4 1.6 1 



(b) 




FIG. 11: Bifurcation diagrams in ADE for cases with 
nonzero cross immunity. For each ADE value, we show max- 
ima (blac k) a nd minima (gray), (a) weak cross immunity. 



0.05. 



(b) strong cross immunity, e — 0.6. (For compari- 



son. Fig. |2j shows the case of no cross immunity.) 



of efficient disease control strategies. 

When considering strong cross immunity, most of the 
parameter regions predict unstable steady state behav- 
ior, as shown in Fig.[9l In fact, when the cross immunity 
parameter e is greater than 0.65, stable endemic behavior 
was achieved only for unrealistically large values of ADE. 
As a result, strong cross immunity destabilized the sys- 
tem, and we observed complicated aperiodic fluctuations, 
such as quasiperiodic behavior and chaotic outbreaks. In 
contrast to the synchronized periodic behavior seen for 
weak cross immunity, we observed that both quasiperi- 
odic and chaotic attractors exhibited strains that were 
unsynchronized. Asynchrony in chaotic outbreaks has 
also been observed in multistrain models with ADE and 
no cross immunity (28j . 

Because time series data for dengue fever show asyn- 
chronous outbreaks for the different strains and non- 
periodic behavior (25l | , our work suggests possible refined 
parameter ranges for dengue in terms of ADE and cross 
immunity. Specifically, either the ADE or the temporary 
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cross immunity must be strong enough to put the sys- 
tem in the chaotic, desynchronized region, where certain 
types of unstable steady states were observed. There is 
now a need to quantify multistrain models against ex- 
isting data sets (such as [25]) and further refine param- 
eter estimates. It should be noted that the model pre- 
sented here does not include seasonality. Because dengue 
is carried by mosquitoes and displays outbreaks with a 
seasonal component, including annual variations in the 
contact rate will likely be necessary for good quanti- 
tative agreement between models and data. However, 
other longer period components exist in the data, and 
are probably due to the interaction between the seasonal 
contact rate fluctuations and the instabilities induced by 
the ADE and cross imunity parameters. From the ADE 
model analyzed in 0, , it was observed that the mean 
period of oscillations was very sensitive with respect to 
the ADE parameter. In the current work, we have also 
done a preliminary sensitivity analysis of the mean os- 
cillation period on the cross immune response parame- 



ter. Here we found that small changes in e may yield 
very large changes in the oscillation period. Therefore, 
in order to connect the model with measured mean pe- 
riods from data, both antibody enhancement and cross 
immunity will play an important role in model predic- 
tion and control. In closing, there are many other mod- 
eling variations which we have omitted but will refine 
model fidelity in future work. These include inhomo- 
geneity in contact rate due to spatial density variation in 
the mosquito populations, fluctuations in the sociologi- 
cal parameters such as contact, birth and death rates, as 
well as general stochastic fluctuations in the population 
itself. Such stochastic effects in finite populations, which 
can lead to fadeout of the disease, may also impact future 
disease controls. 
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